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ABSTRACT 

The precise reconstruction of the turbulent volume is a key point in the development of new-generation Adaptive 
Optics systems. We propose a new profilometry method named CO-SLIDAR (COuplcd Slope and scintillation 
Detection And Ranging) , that uses correlations of slopes and scintillation indexes recorded on a Shack- Hartmann 
from two separated stars. CO-SLIDAR leads to an accurate retrieval for both low and high altitude layers. 
Here, we present an end-to-end simulation of the C„ profile measurement. Two Shack-Hartmann geometries are 
considered. The detection noises are taken into account and a method to subtract the bias is proposed. Results 
are compared to profiles obtained from correlations of slopes only or correlations of scintillation indexes only. 
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1. INTRODUCTION 

An accurate knowledge of the vertical distribution of the strength of atmospheric optical turbulence is a crucial 
point in the development of Adaptive Optics (AO) facilities. New Wide Field Adaptive Optics (WFAO) concepts 
are investigated for the Extremely Large Telescopes (ELTs) in order to increase the field of correction. The impact 
of the C% profile on the performances of WFAO systems has been pointed outP leading to the need for precise 
tomographic reconstruction of the turbulence volume. 

CO-SLIDAfP is a new method to measure a high-resolution profile. The principle is to use correlations 
of wavefront slopes and scintillation indexes recorded with a Shack-Hartmann (SH) from a binary star. CO- 
SLIDAR combines sensitivity to ground and low altitude layers, with correlations of slopes, like a SLODARPbut 
it also provides sensitivity to weak high altitude layers, taking advantage of correlations of scintillation indexes, 
like SCIDARPor MASSP 

CO-SLIDAR first tests in simulatiorP and on-real datcP are encouraging. The extension of the method to 
a single source has also been tested on infrared dataP in an endo-atmospheric context. We are now looking 
towards a complete on-sky validation of the concept. To prepare it, we perform an end-to-end simulation of 
a C\ profile measurement with CO-SLIDAR, in a concrete astronomical case, with binary stars observed on a 
1.5-meter telescope, taking into account realistic fluxes and detection noises. 

In Section [5J we recall CO-SLIDAR theoretical background. In Section [31 we set out the simulation of 
astronomical turbulent images with two SH geometries. Section 2] is dedicated to the data processing. In 
Section [5J reconstruction results are shown and commented. Our conclusions and perspectives are exposed 
in Section [5] 

Further author information: juliette.voyez@onera.fr 



2. PROBLEM STATEMENT 



The analytical formulation and notations are recalled here, for the understanding of the problem statement and 
its inversion, to get an accurate measurement of the C 2 profile with Shack- Hart mann data. We consider two 
stars with separation 9 in the field of view (FOV). The SH delivers a set of wavefront slopes and intensities 
per frame and per star. For one star, we denote s m (a) the slope measured in the subaperture m, a being the 
position of the star. It is a bidimensionnal vector with two components sj,(a), along the k axis (k G {x,y}), 
corresponding to the two "tip" and "tilt" directions. The star intensity is denoted i m (a), and leads to the 
scintillation index 5i m (a) — where o m (a) is the time-averaged star intensity. Slope correlations 

(s^ n s l n )(6) and scintillation index correlations (Si m 5i n )(8) can be written as C^(h) integrals weighted by functions 
denoted W*f and Wu in the following expressions: 

p+oo 

{s k m s l n ){6)= C 2 n (h)W^(Kd mn ,6)dh, (1) 
Jo 

p + QO 

(Si m Si n )(9) = / Cl (h) W u (h, d mn , 9) dh, (2) 
Jo 

These expressions are derived from the terms of the anisoplanatism error under Rytov approximation.^ The 
weighting functions W depend on SH geometry, statistical properties of the turbulence, star separation 9, dis- 
tance between subapertures d m „ and altitude h. They can be seen as the response of the system to a single layer 
at altitude h, for a certain distance between subapertures and a certain star separation 9. Cross-correlations 
combine two directions of analysis corresponding to the binary star separation 9, while auto-correlations cor- 
respond to the case 9 = 0. Slope correlations arc mainly sensitive to ground and low altitude layers whereas 
scintillation correlations are more sensitive to high altitude layers^ but there is no scintillation on the pupil. 
Using cross-correlations, the altitude resolution dh and the maximum sensing altitude H max are obtained with 
simple geometrical rulea^ involving the subaperture diameter d su b and telescope diameter D: 

dh ~ % (3) 
9 

Hmax — ~q (4) 

Experimentally, correlations are estimated from a finite number of frames. Then, they are arranged in a 
single dimension covariance vector C mes , which is directly related to the profile according to Eq. [5] 

C mes = MCl + C d + u (5) 

where M is the interaction matrix with column vectors formed by the concatenation of the weighting functions 
W. Cd is the covariance vector of detection noises affecting slope and intensity measurements, u is the conver- 
gence noise representing uncertainties on C mes due to the limited number of frames. C mes is estimated from 
measurements affected by photon and detector noises that bias the correlation estimates. Assuming the system 
is well calibrated, it is possible to completely determine Cd and define a non-biased estimation of the correlation 
vector, C mes = C mes — Cd- Finally, this makes it possible to rewrite the problem statement: 

C mes = MCl + u (6) 

The covariance matrix of u, C conv = (uu T ), is estimated from C mes and Cd as the empirical covariance matrix 
of a Gaussian random variable vector. A sampled estimate of C^, S, can be retrieved from the inversion of Eq. [5] 
Under positivity constraint, because C 2 is never negative, S minimizes the maximum likelihood (ML) criterion 
J: 

J = (C mes - MSfC-^ (C mes - MS) (7) 

The diagonal of the covariance matrix (M T C~ l lv M)~ 1 can be used as an upper bound of the (square of the) 
sought error bars on the profile S. The implementation of these error bars will cover another material. 



3. SIMULATION OF ASTRONOMICAL TURBULENT IMAGES 



This section is dedicated to the description of an end-to-end simulation in order to produce turbulent SH images 
of a binary star, taking into account photon and detector noises. Parameters of the simulation are given in 
Subsection 13.11 The numerical modeling is described in Subsection 13.21 Simulated images are finally presented 
in Subsection 13.31 

3.1 Simulation parameters 

We consider a D = 1.5 meter telescope with 30 % of central obscuration, and two SH geometries, 30 x 30 and 
15 x 15 subapertures. The subapcrturc diameter is d su b = 5 cm or d su b = 10 cm, depending on the number of 
subapertures. The wavelength is A = 0.55 fim. The observed object is a binary star with separation 9 = 20 " 
modeled by a two-point source. We assume a difference of one magnitude between the two stars. Fluxes are about 
120 and 300 photons per subaperture and per frame for each star. The profile is typical of an astronomical 
site, with strong turbulence at ground level, and turbulent activity in altitude between 13 and 15 km (Fig. [TJ. 
It is composed of 32 C\ values defined at 32 altitudes (i.e. 32 layers). The resultant Fried parameter ro ~ 5 cm. 
Outer scale Lq = 8 m and inner scale Iq = 5 mm. Simulation outputs are 660 x 660 Shannon sampled images, 
leading respectively to 22 x 22 and 44 x 44 pixels per subaperture, for each SH geometry 
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Figure 1. Theoretical C„ profile used for the simulation. 

3.2 Numerical modeling 

The numerical modeling is composed of two parts, the first consists of the diffractive propagation and the 
second contains the image formation process. Here we expose briefly how it works, more details can be found in 
RefP In order to simulate the propagation through turbulence, we use the PILOT (Propagation and Imaging, 
Laser and Optics through Turbulence) code, developed at Onera. The source is sampled on a bidimensionnal 
Cartesian grid of points. The turbulent volume is broken down into a series of discrete layers. The code simulates 
turbulent phase screens, representing turbulent layers, following the von Karman spectrum of the refractive- index 
fluctuations, and Fresnel propagation takes place between the screens. Therefore we consider both phase and 
scintillation effects. The output of PILOT is a complex electromagnetic field at the telescope pupil. It is then 
sampled at SH subaperture level. Point-spread functions (PSF) for each point of the source are then obtained by 
computing the square modulus of the Fourier transform of the electromagnetic field. The final turbulent image 
in a SH subaperture is built with the sum of the PSFs. 

3.3 Shack- Hart mann images 

30 x 30 and 15 x 15 SH images are computed after 100 wave propagations through 32 phase screens, representing 
the 32 turbulent layers of the profile in Fig. [1] Then, 9 cuts are performed on the electromagnetic field after 



each propagation, in order to get 900 SH frames. These images are noise-free. Detection noises are added outside 
of the simulation code. Photon noise is modeled by a Poisson law and detector noise, assumed to be a e - = 1 e" 
per pixel and per frame, is modeled by a Gaussian law. Here, we assume that noise is statistically independent 
for different subapertures and stars. Typical SH frames and subaperture images are presented in Fig. [2j 
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Figure 2. Full SH images (averaged over 100 frames) and short-exposure subaperture images, obtained after propagation 
through turbulence and addition of detection noises. Left: 30 x 30 geometry. Right: 15 x 15 geometry. 

4. DATA PROCESSING 

In this section we describe how we measure slopes and scintillation indexes from simulated SH images in 
Subsection 23] in order to make correlation maps. The bias subtraction process is detailed in Subsection 14.21 

4.1 Slopes and scintillation indexes extraction 

Wavefront slopes and scintillation indexes for each star are calculated in each subaperture for all SH short- 
exposure images. Subapertures that are less than 95 per cent illuminated are excluded from the analysis. To 
compute slopes s m (a) in each of the two orthogonal "tip" and "tilt" directions, we use a centre of gravity (COG) 
algorithm. In order to separate the contributions of the two stars, and to limit the effects of noise, centroids are 
evaluated by delimiting windows around the stars. Boxes of 4 x 4, 6 x 6 and 8x8 pixels are considered for the 
30 x 30 geometry, where the spots are not so much distorted because of turbulence {^f^- ~ 1). For the 15 x 15 
geometry, spots are more distorted. Indeed, the subaperture diameter is two times larger (-^ ~ 2), and images 
are more affected by turbulence. So, we have to chose larger windows, and size of 8x8, 10x10, 12x12, 14x14 
and 16 x 16 pixels are considered. The positions of the boxes arc determined by locating the mean position of 
the maximum of the spots in the subapertures. The total intensity corresponds to the total of all pixel intensities 
into the box. Scintillation index 5i m (a) is derived by subtracting the mean intensity in the subaperture over the 
whole sequence of images o m (a), and by dividing the result by this same term. Slopes and scintillation indexes 
are extracted from noise-free and noisy images for comparison. Noise propagation occurred in the computation 
of slopes and scintillation indexes. As they are measured using a COG algorithm, noises contribution are then 
given byP 

TT 2 1 f N T \ 2 



for slope measurements. Nt is the image full width at half maximum (FWHM), Njj is the FWHM limited by 
diffraction and N$ is the number of pixels for the COG calculation. However, in Eq. [SJ t does not depend 

on the size of the window. It has been shown that when TVs > 2JVu, h t can be written asP^ 

N ph \N D 



For scintillation measurements we have: 



1 



(11) 



ph 

For each kind of data, the total noise, &Ad, naiae for slopes and o'gi noise for scintillation, is the sum of the terms <r^ hot 
and cr^ et . The effects of the size of the window and the effects of noise on the measurement are investigated 
in Subsections 15.11 and 15.21 

4.2 Making of correlation maps and bias subtraction process 

Auto-correlations (on a single star) and cross-correlations (between the two stars) are computed from slopes and 
scintillation indexes. In this simulated case, auto-correlations are evaluated for each star, and averaged, which 
increases the statistics by doubling the number of samples. Correlations arc calculated for all separations between 
subapertures and represented as correlation maps. These maps show the mean correlation between all couples 
of subapertures having the same gap (Fig. [3]), for the two kind of data, slopes and scintillation indexes. In the 
presence of measurement noises, these maps are biased. Assuming that the noise is statistically independent 
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Figure 3. Correlation maps for the 30 x 30 geometry. The maps dimensions are 2 x number of subapertures- 1, illustrating 
the correlation for all couples of subapertures. Top-left: slope auto-correlations in the x direction. Top-right: slope cross- 
correlations in the x direction. Bottom-left: scintillation auto-correlations. Top-right: scintillation cross-correlations. 



for different stars and subapertures, only the variance, namely the central value of the auto-correlation maps of 
slopes and scintillation are biased. A good way to estimate the bias is then to analytically calculate the noise 
using Eqs. [5J [HI HH1 for slope auto-correlations and Eqs. [11] [TJJ for scintillation auto-correlations, knowing 
the number of photons N p h and the detector noise a e - . As the two stars have different magnitudes, we apply 
different numbers of photons per subaperture and per frame in the formulae. The final bias is the mean of 
the biases given by each star. Once it is calculated, we just have to subtract it from the central value of the 
auto-correlation maps. The effect of the bias subtraction on the retrieval will be studied in Subsection 15.21 

5. C 2 N RESTORATION: RESULTS 

The C% profiles reconstructed with the two SH geometries are presented and discussed in this section. We begin 
by studying the influence of the size of the window in Subsection 15.11 Then, the effects of measurement noises 
and bias subtraction are analyzed in Subsection 15.21 We finally compare the results with those obtained with 
other methods in Subsection 15.31 and the advantages of the two SH geometries are examined. 

5.1 Influence of the COG window 

As we said in Subsection 14.11 windows of different sizes are considered to compute slopes and intensities. Here, 
we only consider noise-free data. C\ profiles reconstructed with the two geometries, for different sizes of window, 
are shown in Fig. [4] values lower than 1 x 10 -19 m~3 are automatically put to 1.3 x 10~ 19 m~s on the 
graph, for better understanding of the estimation. With the 30 x 30 geometry we restore 32 layers, that is to 
say the nominal number of layers of the input theoretical profile, while with the 15 x 15 geometry, we only 
restore 20 layers. We can notice immediately that the size of the window introduces a slight bias in the profile 
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Figure 4. Impact of the size of the window for COG and intensity calculation on the reconstructed C„ profile. Left: 30 x 30 
geometry. Right: 15 x 15 geometry. 

measurement, and that for the two SH geometries. For the 30 x 30 geometry, although the low altitude layers 
are almost well reconstructed with all sizes of window, the greatest differences occurred after 4 km of altitude. 

values are always over-estimated, but the smaller is the box, the higher is the over-estimation. The layer at 
14 km is imperfectly reconstructed with the 6x6 and 8x8 windows, while it is not at all with the 4x4 window. 
Layers after 15 km of altitude are not well estimated with the three sizes of window (c/. Subsection 15. 3p . For 
the 15 x 15 geometry, we also observe that the reconstruction is improved when the size of the box increases. 
Identically, low altitude layers are almost well reconstructed. But, the plateau between 6 km and 12 km is better 
estimated when the box dimensions reach 10 x 10 pixels. Results are quite similar for all sizes of window for 
layers between 12 km and 14 km. The layer at 15 km is better evaluated when the box size reaches 12 x 12 
pixels. As for the other geometry, layers after 15 km of altitude are not well estimated with all sizes of window. 
Windows of 4 x 4 and 6x6 pixels were also tested, but results were too much biased so they are not presented. 

This bias depending on the size of the window can be explained. When calculating slopes and intensities 
in little boxes, smaller than the whole subaperture, the spot is partially truncated. So, part of the flux is not 



taken into account for the COG calculation and the total intensity evaluation, leading to a biased estimation, 
that propagates in the C\ profile retrieval. As a consequence, the size of the window must be chosen carefully. 
It must depend on the spot distortion and the signal to noise ratio (SNR) in the subaperture. The more the spot 
is distorted, the more the box has to be large, in order not to truncate too much the spot for the measurements. 
However, a larger window introduces more noise in the measurements (Eqs. l9l fTOl IT2|) . Moreover, if stars are too 
closed, the size of the box is implicitly limited. In the following, we keep windows of 8 x 8 pixels for the 30 x 30 
geometry and windows of 12 x 12 pixels for the 15 x 15 geometry. 

5.2 Influence of noise and effect of bias subtraction 

Now we perform reconstruction from noisy data. The impact of noise propagation on the results is presented 
in Fig. [5l Noise adds an offset to all estimated values, with both geometries. Low altitude layers, until 4 km 




of altitude, are less impacted by this offset than layers between 4 and 15 km. Actually, these layers are mostly 
restored thanks to the scintillation signal. But, this signal has a poorer SNR than the slope signal. For the two 
geometries, the SNR on slopes and scintillation signals are detailed in Tables [T] and [21 using Eqs. M ITUI UTI [T21) . 
These tables show that there is a good SNR on slope signal <j\^ t b , but it is very weak on the scintillation signal 
a Si t b > even smaller than 1 in the case of the 15 x 15 geometry. This observation explains the bad estimation of 
the profile at altitudes between 4 and 15 km. 
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Table 1. SNR on slope signal, for the two SH geometries. 
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Table 2. SNR on scintillation signal, for the two SH geometries. 



Then, we carry out the bias subtraction, as explained in Subsection 14.21 and using and cr^ 

values presented in Table [1] and They arc subtracted from the central value of the auto-correlation maps. 
Bias-subtracted correlations are then used to restore the profile, and results are presented in Fig. El This 




correction nearly allows to restore the unbiased profile, especially with the 30 x 30 geometry. Slight differences 
can be noticed with the 15 x 1 5 geometry, probably due to the poor SNR, but the efficiency of the method is 
demonstrated. 

5.3 Comparison with other methods, advantages of the two SH geometries 

We finally compare profiles restored with CO-SLIDAR with C\ profiles restored with slope data only or 
scintillation data only. Results are shown in Fig. [7] Comments are common to both geometries. In all cases, 
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Figure 7. Comparison of CO-SLIDAR with other methods. The bias has been subtracted in the three cases. Left: 30 x 30 
geometry. Right: 15 x 15 geometry 



noisy data were used, and the bias has been calculated and subtracted. Slope data allow a good reconstruction 
of ground and low altitude layers, but they provide worse sensitivity at high altitude. Scintillation data perform 
a good reconstruction at low and high altitude, but do not permit to estimate the ground layer, as there is no 
scintillation on the pupil. This goes in the way of using both type of data to retrieve the C% profile. Actually, 
CO-SLIDAR combined the advantages of the two kind of data in a single instrument, and allows to estimate 
an accurate C\ profile from the ground to 15 km of altitude, with these SH geometries. Layers after 15 km 
of altitude are not well reconstructed, neither with CO-SLIDAR nor slope nor scintillation data only. Indeed, 
after 15 km, with these geometries, cross-correlations are blind to turbulence, according to Eq. SJ Information 
can then be provided by scintillation auto-correlations, but here the scintillation signal seems to be too weak to 
perform a good estimation of values. Better sensitivity could be achieved using smaller subaperturesj^ but 
this would lead to deal with very low fluxes at subaperture level in an astronomical context. Anyway, these 



values are very small, < 1 x 10~ 18 m _ S i and represent a negligible part of the whole turbulence. Nevertheless, 
these layers over 15 km of altitude could be estimated using a binary star with a smaller separation, to increase 
the altitude range sensitivity, but this would decrease the altitude resolution. 

Finally, the performances of the two SH geometries need to be compared. Reconstructed profiles are presented 
in Fig. [HI We immediately notice the 30 x 30 geometry has a better altitude resolution, about 500 m, while the 
15 x 15 geometry only permits a resolution of about 1 km, as forecast by Eq. [31 Both SH geometries allows a 
profile measurement from the ground to 15 km but do not permit to go over with reliable results, as explained 
before. The 15 x 15 geometry, by collecting more flux, allows a better sky-coverage. Nevertheless, one has to be 
careful because if turbulence is too strong, images will be formed of speckles, leading to a difficult extraction of 
slopes and scintillation indexes. Both geometries have their advantages and they could be used together in order 
to cover a full observation night. 
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Figure 8. Comparison of the two SH geometries. 



6. CONCLUSIONS AND PERSPECTIVES 

In this paper, we detailed the performances of CO-SLIDAR in an end-to-end simulation, in a realistic astronomical 
case. Detection noises have been taken into account, and a way to subtract the bias from the measurements has 
been proposed and its efficiency demonstrated. CO-SLIDAR has been compared to other methods, highlighting 
the good complementarity of correlations of slopes and scintillation indexes, to provide turbulence sensitivity at 
low and high altitude. Two SH geometries have been proposed, each one with its proper advantages. 

In order to perform complete on-sky validation, we acquired CO-SLIDAR data through an observation cam- 
paign on the MeO 1.5-meter telescope, with the two SH geometries, at the Cote d'Azur Observatory, in South 
of France. Data processing is in progress, in order to estimate high-resolution C\ profiles. SCO-SLIDAR, the 
extension of the method to a single source, should also be tested on single star data. Computation of error bars 
on the reconstructed profile has to be implemented, together with a refined study of the convergence noise. As 
sensitivity to the outer scale Lq has been proved in a previous workpl combined estimation of the two parameters 
is conceivable. 
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